MODULE fefmm_operator

   USE fmm_global_paras
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: fefmm_build_operator, fefmm_get_FEop

   REAL(REALK), ALLOCATABLE :: FEop(:,:)

   REAL(REALK), PARAMETER :: T_VALUE_THR = 0.5D0

   INTEGER(INTK) :: LIPN
   INTEGER(INTK) :: FEdim, FEdimB, NFS
   REAL(REALK)   :: FEgrain

CONTAINS

!-------------------------------------------------------------------------------

   SUBROUTINE fefmm_build_operator(scheme)

      IMPLICIT NONE
      TYPE(scheme_paras), INTENT(IN) :: scheme


      REAL(REALK), PARAMETER :: x(5) = (/0.1488743389d0, 0.4333953941d0, &
                                         0.6794095682d0, 0.8650633666d0, &
                                         0.9739065285d0/)
      REAL(REALK), PARAMETER :: wt(5) = (/0.2955242247d0, 0.2692667193d0, &
                                          0.2190863625d0, 0.1494513491d0, &
                                          0.0666713443d0/)

      REAL(REALK) :: tmax = 10000.0d0
      REAL(REALK) :: ncell = 25

      REAL(REALK)   :: FEgrain
      REAL(REALK)   :: t1,t2
      REAL(REALK)   :: xr,xm,dx
      INTEGER(INTK) :: tn, np,nc

      print *, "Building Coulomb operator..."

      FEdim   = scheme%FEdim   ! number of grid points in one box
      FEdimB  = (3*FEdim -2)   ! number grid points in near-field space
      NFS     = 3              ! number boxes in near-field space

      LIPN    = scheme%lipn
      FEgrain = scheme%grain / (FEdim -1)

      np = (FEdimB)**3    ! rho-point
      nc = FEdim**3          ! V-point
!      nc = 1
      ALLOCATE (FEop(np,nc))
      print '(A,F10.2,A)', " Memory for FEop = ", SIZE(FEOp)*8d0*1d-6, " MB"

      FEop(:,:) = zero
! Gauss-Legendre 10-point quadrature over hierarchical cells
      DO nc = 1, ncell
         t2 = tmax*(0.5d0**(nc-1))
         t1 = t2*0.5d0
         IF (nc == ncell) t1 = 0d0
         xm = 0.5d0*(t2+t1)
         xr = 0.5d0*(t2-t1)
         DO np = 1, 5
            dx = xr*x(np)
            CALL fefmm_update_FEop(xm+dx,wt(np)*xr)
            CALL fefmm_update_FEop(xm-dx,wt(np)*xr)
         END DO
      END DO

      FEop(:,:) = FEop(:,:)*two/ROOTPI

      print *, "Storing Coulomb operator..."
      CALL fefmm_store_FEop

      DEALLOCATE (FEop)

   END SUBROUTINE fefmm_build_operator

!-------------------------------------------------------------------------------

   SUBROUTINE fefmm_update_FEop(t,fac)

      USE fefmm_LIP1_Fij, ONLY: fefmm_build_Fij_LIP1
      USE fefmm_LIP2_Fij, ONLY: fefmm_build_Fij_LIP2
      USE fefmm_LIP3_Fij, ONLY: fefmm_build_Fij_LIP3

      IMPLICIT NONE
      REAL(REALK), INTENT(IN) :: t, fac

      REAL(REALK) :: Fij(FEdimB,FEdim)

!print *, 'lambda, t, fac =', FEgrain, t, fac
      SELECT CASE (LIPN)
         CASE (1)
            CALL fefmm_build_Fij_LIP1(FEdim,NFS,FEgrain,t,Fij)
         CASE (2)
            CALL fefmm_build_Fij_LIP2(FEdim,NFS,FEgrain,t,Fij)
         CASE (3)
            CALL fefmm_build_Fij_LIP3(FEdim,NFS,FEgrain,t,Fij)
         CASE DEFAULT
            CALL fmm_quit('LIPN too large in build FEop!')
      END SELECT

      CALL fefmm_contract_Fxyz(fac,Fij)

   END SUBROUTINE fefmm_update_FEop

!-------------------------------------------------------------------------------

   SUBROUTINE fefmm_contract_Fxyz(fac,Fij)

      IMPLICIT NONE
      REAL(REALK), INTENT(IN) :: fac, Fij(:,:)

      INTEGER(INTK) :: i,j,k, u,v,w, ijk, uvw

      ijk = 0
      DO i = 1, FEdim        ! V-point
      DO j = 1, FEdim
      DO k = 1, FEdim
         ijk = ijk +1
         uvw = 0
         DO u = 1, FEdimB    ! rho-point
         DO v = 1, FEdimB
         DO w = 1, FEdimB
            uvw = uvw +1

            FEop(uvw,ijk) = FEop(uvw,ijk) + fac*Fij(u,i)*Fij(v,j)*Fij(w,k)

         END DO
         END DO
         END DO

      END DO
      END DO
      END DO

   END SUBROUTINE fefmm_contract_Fxyz

!-------------------------------------------------------------------------------

   SUBROUTINE fefmm_store_FEop

      IMPLICIT NONE

!fixme
      CHARACTER(LEN=11), PARAMETER :: Name = 'FE_operator'
      CHARACTER(LEN=255) :: FBuf
      INTEGER(INTK) :: i,j,k, u,v,w, ijk, uvw

      FBuf = TRIM(Name)//".fefmm"
      OPEN(UNIT=IO2, FILE=TRIM(FBuf), STATUS='REPLACE',  &
           ACCESS='SEQUENTIAL', FORM='UNFORMATTED')
      REWIND(IO2)

      ijk = 0
      DO i = 1, FEdim        ! V-point
      DO j = 1, FEdim
      DO k = 1, FEdim
         ijk = ijk +1
         uvw = 0
         DO u = 1, FEdimB    ! rho-point
         DO v = 1, FEdimB
         DO w = 1, FEdimB
            uvw = uvw +1

            WRITE(IO2) FEop(uvw,ijk)

         END DO
         END DO
         END DO
      END DO
      END DO
      END DO

      CLOSE(UNIT=IO2, STATUS='KEEP')

   END SUBROUTINE fefmm_store_FEop

!-------------------------------------------------------------------------------

   SUBROUTINE fefmm_get_FEop(FEop)

      IMPLICIT NONE
      REAL(REALK), INTENT(OUT) :: FEop(:,:)
!fixme
      CHARACTER(LEN=11), PARAMETER :: Name = 'FE_operator'
      CHARACTER(LEN=255) :: FBuf
      INTEGER(INTK) :: i,j,k, u,v,w, ijk, uvw

      FBuf = TRIM(Name)//".fefmm"
      OPEN(UNIT=IO2, FILE=TRIM(FBuf), STATUS='OLD', ACTION='READ',   &
           ACCESS='SEQUENTIAL', FORM='UNFORMATTED')
      REWIND(IO2)

      ijk = 0
      DO i = 1, FEdim           ! V-point
      DO j = 1, FEdim
      DO k = 1, FEdim
         ijk = ijk +1
!         ijk = 1
         uvw = 0
         DO u = 1, FEdimB    ! rho-point
         DO v = 1, FEdimB
         DO w = 1, FEdimB
            uvw = uvw +1

            READ(IO2) FEop(uvw,ijk)
!            print *, uvw, ijk, FEop(uvw,ijk)

         END DO
         END DO
         END DO
      END DO
      END DO
      END DO

      CLOSE(UNIT=IO2, STATUS='KEEP')

   END SUBROUTINE fefmm_get_FEop

!-------------------------------------------------------------------------------

END MODULE fefmm_operator
